### Load Data
dat <- readRDS("./dat.rds")
str(dat)
head(dat)
dat[,1:55] = dat[,1:55]
dat[,57] = dat[,57] 
### Function to calculate error

rsquared <- function(y,f)
{
  return (1- sum( (y-mean(y))^2 )/sum( (y-f)^2 ) )
}

mase <- function(y,f)
{
  return(mean(abs( (y-f) / mean(abs(y-mean(y))))))
}

mae <- function(y,f)
{
  return( mean(abs(y-f)/y) )
}
###------
### Correlation analysis
library(caret)
dat.cor <- cor(dat[,1:55])
summary(dat.cor[upper.tri(dat.cor)])
highcor <- caret::findCorrelation(dat.cor,cutoff=0.9)

###------
### Bootstrapping
bsize = 100
sampleSize = 25*100
responseLocation = 56
trainSize = 0.7
betaMatrix <- matrix(ncol=bsize,nrow=56)
mseMatrix.train <- matrix(ncol=bsize,nrow=1)
mseMatrix.test <- matrix(ncol=bsize,nrow=1)

for (i in 1:bsize)
{
  trainObs <- floor(trainSize*dim(dat)[1])
  testObs <- trainObs + 1
  idx <- sample(c(1:33),size=33,replace=F)
  trainIDX <- idx[1:trainObs]
  testingIDX <- idx[trainObs:dim(dat)[1]]
  
  training <- dat[trainIDX,]
  x.test <- as.matrix(dat[testingIDX,1:55])
  y.test <- as.matrix(dat[testingIDX,responseLocation])
    
  bs <- sample(c(1:trainObs), size=sampleSize,replace=T)
  datMod <- training[bs,]
  dim(datMod)
  
  ###------
  
  x.train <- as.matrix(datMod[,1:55])
  y.train <- as.matrix(datMod[,responseLocation])

  x.train <- scale(x.train)
  x.test <- scale(x.test)
  datCaret <- datMod[,c(1:55,responseLocation)]
  
  library(glmnet)
  mod1 <- cv.glmnet(x=x.train[,1:55], y= y.train, family="gaussian",alpha=1, nlambda=100)
  
  mod2 <- glmnet(x=x.train[,1:55], y= y.train, family="gaussian", alpha=1, lambda=mod1$lambda.1se)
  beta <- as.matrix(coefficients(mod2))
  betaMatrix[,i] <- beta
  
  
  yhat.test <- predict(mod2, s=mod1$lambda.1se, newx=x.test)
  yhat.train <- predict(mod2,s=mod1$lambda.1se, newx=x.train)
  mseMatrix.train[1,i] <- mae(y.train,yhat.train)
  mseMatrix.test[1,i] <- mae(y.test,yhat.test)
}

plot(mseMatrix.test[1,])
abline(h=mean(mseMatrix.test),col='red')

betaCount <- ifelse(betaMatrix==0,0,1)
importance <- as.matrix(apply(betaCount,1,sum))

library(ggplot2)
library(magrittr)
slct <- data.frame(Variable = paste("X",c(1:55),sep=''),Freq=importance[2:56,1])
slct %>%
  ggplot() +
  geom_bar(aes(x=Variable, y=Freq), stat='identity', width=0.1) +
  geom_point(aes(x=Variable, y=Freq)) +
  coord_flip()

slctVariable <- data.frame(Variable=slct[slct$Freq>1,1], Location=which(slct$Freq>1)) 
slctVariable

mseDF <- data.frame(Train=t(mseMatrix.train),Test=t(mseMatrix.test))
boxplot(mseDF)

saveRDS(slctVariable,file="./slctVariable.rds")
save.image("./FactorImportance.RData")
